The Type II Phase Resetting Curve is Optimal for Stochastic Synchrony 
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The phase-resetting curve (PRC) describes the response of a neural oscillator to small perturba- 
tions in membrane potential. Its usefulness for predicting the dynamics of weakly coupled determin- 
istic networks has been well characterized. However, the inputs to real neurons may often be more 
accurately described as barrages of synaptic noise. Effective connectivity between cells may thus 
arise in the form of correlations between the noisy input streams. We use constrained optimization 
and perturbation methods to prove that PRC shape determines susceptibility to synchrony among 
otherwise uncoupled noise-driven neural oscillators. PRCs can be placed into two general categories: 
Type I PRCs are non-negative while Type II PRCs have a large negative region. Here we show that 
oscillators with Type II PRCs receiving common noisy input sychronize more readily than those 
with Type I PRCs. 



Introduction 

Synchronous oscillations are found in many brain areas 
and are responsible for macroscopic electrical responses 
of the brain including field potentials and EEG signals. 
Within a single brain area, synchronization of neuronal 
activity serves to amplify signals to upstream regions [l[ , 
while synchronization across different areas may allow 
activity to be selectively routed. 

Considerable theoretical interest has recently emerged 
in the generation of synchrony by correlated "noisy" in- 
puts to uncoupled oscillators 0, d, 0, || , a phenomenon 
we will refer to as stochastic synchrony. In the brain, 
stochastic synchrony may account for observations such 
as long-range synchronization 0, 0] , that are difficult to 
explain by the presence of synaptic connectivity alone. 
Moreover, noisy inputs have been shown to synchronize 
real neurons in vitro @. 

The key component in the study of noisy oscillators is 
the phase-resetting curve (PRC). This curve character- 
izes how inputs to an oscillator shift its timing, or phase. 
In the context of neurons, spike times are believed to play 
an important role in coding and in the propagation of in- 
formation across brain regions. Thus, the PRC provides 
a quantitative characterization of how inputs to neural 
oscillators alter the timing of spikes. 

The theory of deterministic oscillators has shown that 
the type of bifurcation from steady-state to periodic be- 
havior determines the shape of the PRC. Weak coupling 
theory shows that the form of the interaction between os- 
cillators together with their intrinsic response (the PRC) 
provide sufficient information about the ability of the 
coupling to synchronize (or desynchronize) the oscilla- 
tions. For very fast excitatory synaptic interactions, 
Type II oscillators characterized by the Hopf bifurcation 
synchronize more readily than Type I oscillators charac- 
terized by the saddle-node-on-an-invariant-circle (SNIC) 
bifurcation [l(J HH, E3 ■ This difference in ability to 
synchronize with excitatory coupling is a consequence of 
the shape of the PRC occurring near the two different 
bifurcations. A PRC which contains both negative and 
positive lobes can allow inputs to both slow down the os- 



cillator which is ahead and speed up the oscillator which 
is behind. In contrast, a non-negative PRC can only 
speed up the timing of both oscillators, so that synchro- 
nization becomes more difficult. A number of authors 
EEHIl have shown that the PRC near a SNIC is 
non-negative and approximately proportional to 1 — cosi, 
while the PRC near a Hopf is proportional to sin(i + a). 
Thus, Type II PRCs have a large negative lobe, whereas 
Type I PRCs are strictly positive. 

Two recent papers have shown that Type II PRCs are 
better than Type I PRCs at synchronizing uncoupled os- 
cillators with correlated input [la, |l6[. That is, for a 
given input correlation of the noisy stimulus, the output 
correlation of the oscillators is higher with Type II than 
with Type I PRCs. In these two papers, specific functions 
for PRCs were checked (namely, sin(t) and 1 — cos(t)), 
and the correlations and degree of synchrony were ana- 
lytically and numerically computed. However, it is not 
known whether there are other PRC shapes that might 
produce even stronger stochastic synchronization. 

The easiest way to quantify stochastic synchrony is 
to examine the Lyapunov exponent, the rate at which 
two oscillators receiving identical inputs converge to syn- 
chrony. In this paper we will explore how this quantity 
depends on the shape of the PRC. In particular, we find 
that Type II PRCs lead to faster convergence than do 
Type I, and we use variational principles to determinine 
the optimal shape of the PRC to maximize this conver- 
gence. 

First in Section [J we introduce the phase reduction 
of a stochastically driven neural oscillator using the Ito 
change of variables, and in Section [IT] we derive the Lya- 
punov exponent for two such oscillators receiving com- 
mon noise. Next we use the Fokker-Planck equation 
in Section IIIII to obtain the probability distribution of 
the phase of a noise-driven neural oscillator. The Euler- 
Lagrange method for constrained optimization allows us 
in Section IIVI to find the PRC that minimizes the Lya- 
punov exponent. This leads to a 4th order system of non- 
linear differential equations, which we approximate to an 
arbitrary order of accuracy using regular perturbations 
in Section [V] The resulting approximation shows that a 
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Type II PRC achieves the minimal Lyapunov exponent, 
hence producing more robust convergence to synchrony 
than a Type I PRC. Several interesting cases that arise 
as a function of the constraint parameters are discussed 
in Section IVII Finally in Section IVIII we show that nu- 
merical solution of the 4th order system agrees with the 
perturbation-derived approximation. 



under some reasonable circumstances the additional term 
can be made arbitrarily small. Thus we will stay with 
the conventional phase-reduced model as first proposed 
by Teramae and Tanaka 

II. LYAPUNOV EXPONENT 



I. 



ITO PHASE REDUCTION 



Consider a neural oscillator with additive white noise 
decribed by the stochastic differential equation 



dX = F(X)dt + aMdW, 



(1) 



where F{X) represents the deterministic equations of 
motion, er is the amplitude of the noise, M is a con- 
stant matrix, and dW is a vector of Gaussian white noise. 
Note that for a general limit-cycle oscillator, there need 
be no constraints on the entries of M. For neural mod- 
els however, the noise typically occurs in current felt by 
the neuron, and this current appears only in the voltage- 
component of the deterministic model. Without loss of 
generality, we take the voltage to be the first component. 
Thus, we will assume here that M has all zero entries 
except for the (1,1) element, which is identically 1. 

The phase reduction method 0] applied to Eq. (TTJ) gives 
a stochastic differential equation for the evolution of the 
oscillator's phase: 



d6 = dt + aA(9)dW, 



(2) 



where we have assumed without loss of generality that 
the intrinsic frequency of the oscillator is lu — 1, and 
dW is now a scalar white noise process. Here A is the 
infinitesimal phase response curve defined by 



A(0) := V x 9 



X=X (6) 



where Xq(9) is the unperturbed limit-cycle solution of 
the deterministic equation X = F(X). Sec Kuramoto 
[13, pages 26-27. 

It is now important to note that the usual phase reduc- 
tion method uses the conventional change of variables, 
so Eq. ^ must be regarded as a Stratonovich differential 
equation 0, [l8| . To eliminate the correlation between 9 
and the white noise £ = dW, we must apply Ito's Lemma 
to obtain an equivalent but analytically more convenient 
formulation 



d6 



-A'(0)A(0) 



dt + aA{6)dW, (3) 



where ' denotes -jL. In a recent paper, Yoshimura and 
Arai [lj| show that Eq.Q is incomplete and that an- 
other term must be added in the case where the noise is 
strictly white. However, more recently (in preparation) 
we show that the correct reduction is more subtle, and 



As a standard measure of susceptibility to synchrony, 
we will now derive the Lyapunov exponent for two identi- 
cal uncoupled neural oscillators receiving common addi- 
tive white noise. The resulting analysis, however, applies 
equally well to an arbitrary number of identical nonin- 
teracting oscillators. 

Let us define the phase difference <j> := 62 — #1, where 
61 and #2 each obey Eq.Q. Linearizing around the syn- 
chronous state (j) — 0, we obtain as in [2J: 



[(A' A)' (0)4>] dt + a[A' 



dW, 



where 9 obeys Eq.Q as well. Since the Lyapunov ex- 
ponent is defined as A := lim^oo l2§i^wj ; j e t U s make 
the change of variables y := \og(<p). Once again we in- 
voke Ito's Lemma, and after simplification we find that 
y satisfies the stochastic differential equation 



dy 



\A"A]dt + aA'dW. 



Next we integrate, divide by t and take the limit as t 
00 to obtain an expression for A. 



A 



t^oo t 

lim £ f A"(9(s))A(6(s))ds 

t^oo It J Q 

ft 



A'(6(s))dW(s) 



Assuming the system is ergodic, we can replace the long 
time average on the right hand side with the spatial or 
ensemble average. Due to the ltd change of variables, the 
last term drops out leaving 



A = y j\"(9)A(6)P(9)d9, 



(4) 



where P{9) is the steady-state distribution of the phase. 

Note that Teramae and Tanaka derive an expression for 
A in Q by making the approximation P{9) = 1. Substi- 
tuting this value into Eq.(j4]) and performing integration 
by parts, they obtain 



j\A'{9)fdB. 



In this paper, however, we wish to retain the generality 
of P{9) as discussed below. 
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III. STEADY-STATE PHASE DISTRIBUTION 



as above and use integration by parts to obtain 



In order to evaluate the Lyapunov exponent, we need 
to obtain the stationary density of the phase when per- 
turbed by noise. Teramae and Tanaka 0] have treated 
the density as uniform, which is correct for weak noise. 
However our subsequent perturbation analysis will re- 
quire higher-order terms, so we will need to derive a more 
accurate value for the steady-state phase distribution. 

By applying the Fokker-Planck equation to , we ob- 
tain after simplification a partial differential equation for 
the probability distribution P(9,t): 



dP 
~dt 



dP 
~d9 



a^8_ 
Yd9 



A 



d(AP) 
39 



Now we may set ^ = to find the steady state, then 
integrate once with respect to 9 to obtain: 



J 



A 



d(AP) 



09 



(5) 



where — J is a constant of integration. We require that 
P(0) = -P(l) and that the solution be normalized, namely 
J P{9) d9 = 1. Note that the equations are singular, 
since A(9) generally vanishes at several places, in par- 
ticular at 9 — 0, 1. In the appendix below, we prove the 
existence of the stationary density by directly solving the 
linear equations and taking appropriate limits. 

In the remainder of this section, we use regular pertur- 
bation theory to approximate the stationary density for 
small noise, < u -C 1. To approximate both J and P 
we substitute 

J = l + a 2 J 1 +<r 4 J 2 + --- 
P{9) = l + <j 2 P 1 (9) + <7 i P 2 (9) + --- 

into equation ([5]). Equating like powers of a gives 



—Ji = -P 1 (9) + -A(9)A'(9). 

Integrating both sides over [0, 1] leaves the constant on 
the left hand side unchanged. For the right hand side, 
note that J* P{6)d6 = 1, and hence J* P 1 (9)d9 = 0. 
Furthermore, A A' = ±^(A 2 ) so that 

Ji = ~(A(1) 2 -A(0) 2 ) 
= 0, 

since A is periodic. Thus we have P x {9) = \A{9)A'{9). 
Similarly, 

-Ji = -P 2 (6) + \a{9) 2 A'{9) 2 + ^A(9fA"(9). 
Since L P2(9)d9 = as well, we can integrate both sides 



■h - 
P2(0) -- 

In summary, 



(A{9)A'(9)) 2 d9 
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^A(60 2 A'(6>) 2 + -A(9) 3 A"(9) 



(A(9)A'(9)) 2 d9. 



J = 1 + ^ j\A{9)A'{9)fd9 



P(9) 



1 + -A{9)A'{9) , 4 



2A{9) 2 A'(9) 2 



-A{9) 3 A"{9)+ / {A{9)A'{9)) 2 d9 
Jo 



(6) 



For the perturbation expansions in the next section, it 
will suffice to write J = 1. We will use Eq.© in Section 
IVII and for the numerical verifications in Section I VIII 



IV. 



CONSTRAINED OPTIMIZATION 



The Euler-Lagrange variational technique provides a 
method for determining the phase resetting curve A that 
minimizes the Lyapunov exponent, subject to appropri- 
ate constraints. To ensure smooth solutions and to elim- 
inate uninformative harmonics of the optimal solution, 
we begin by imposing the general constraint 



a{A(9)) 2 + b{A'{9)) 2 + c(A"(9)) 2 d9 = 1, (7) 



where a, b and c are free parameters. A standard nor- 
malization has a = 1, b = 0, c = 0, but non-zero values of 
b, c endow solutions with additional smoothness. Below 
we will explore the cases that arise from specific choices 
of these. 

We proceed by placing Eqs. ([!]), © and (O together 
with the approximation J = 1 into the Euler-Lagrange 
formula to obtain the functional 

Jo 1 A"AP + v t [aA 2 + b{A') 2 + c(A") 2 - l] 



1 - P + —A(AP)' 



d9 = 0, 



(8) 



where v\ is a free parameter, and v-i (9) represents a con- 
tinuum of free parameters. 
Define the operator 

£(A) := A"AP + Vl [aA 2 + 6(A') 2 + c(A") 2 -l] 



-va(0) 



i -p + —A(Apy 
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The optimal A we seek will satisfy the two equations 



dL d dL d 2 dL 



dA d9dA' d9 2 dA" 
dL d dL 



= 
= 0. 



(9) 
(10) 



dP d6 BP 

Note that we can write two more Euler-Lagrange equa- 



tions, but J^jj- 



simply restates Eq.©, and — 



returns Eq.Q governing P. 



dv 2 



Assuming the parameter c is nonzero, we obtain from 
Eqs.© and (fl"0[) a 4th order system of ordinary differen- 
tial equations: 



P"A + 2{P'A' + PA" + aAv l -bA"v 1 + cA^v 1 ) + ^A{P'v 2 -Pv' 2 )<7 2 = 

AA" - u 2 - ^A(AV 2 + Av' 2 )a 2 = 0. 



(11) 
(12) 



r 



If c = 0, we will have instead the 2nd order system which 
obtains by setting c = in Eq. (fTT|) . When we examine 
the effects of varying the constraint parameters in Section 
VII we will see that the main result remains the same in 
this case as well. 



V. PERTURBATION APPROXIMATION 

Let us first consider the 4th order case where the pa- 
rameter c is nonzero. 

Assuming the noise amplitude a is sufficiently small, 
we write the following expansions 

P (9)+a 2 P 1 (9) + ... 



P(9) 



(13) 



A(0) = A (9)+a 2 A 1 (9) + ... 

V\ = Vlfl + (T 2 V lt l + ... 
V2{0) = ^2,0(#)+^2,lW + - 

Substituting these into (TTTj) and (Tj"2"]) and equating like 
powers of a gives to lowest order: Pq(9) = 1, z^ 2 ,o (0) = 
Aq(9)Aq(6) and the fourth order homogeneous equation 

au h0 A + (1 - 6^, )Ag + c^, A 4) = 0. (14) 
For convenience let us define the differential operator 



d 2 

J = avi fi + (1 - bv l,o)gg2 + CVl ® 



o 4 



Thus Eq.JTJJ becomes J{Aq) — 0, and the first order 
correction Ai obeys the inhomogeneous equation 



J{A 1 



+A (au hl + 3A' A%) 



c^A^. 



Furthermore, substituting the expansions ([13 
Eq-© gives the corresponding constraints: 

aA^ + 6(A^) 2 + c (A^) 2 = 1 

f aA Ai + bA' Q A[ + cAq A'/ = 0. 
Jo 



(15) 
into 

(16) 
(17) 



Before solving Eq. (ri4]) . we must first determine the 
unknown parameter 2/x,o- Since we seek only periodic 
solutions, we can impose a condition on the characteristic 
equation of ifH 



(1 



bvi,o)y 2 + cvi^y A 



0. 



(18) 



Specifically, by requiring that the roots of this polynomial 
satisfy y — 2iri, we determine that 



4^ 2 



^1.0 = 



Abu 2 + 16ctt 4 



Now we are ready to impose periodic boundary con- 
ditions, and we find that the solution of ([14]) is just 
Aq(9) — Co sin(27r(9). The constant of integration Co 
is determined from the constraint (|16[) so that 



Co =± 



V2 



Va + Abu 2 + 16c7r 4 



While both values of Co will give the same minimal value 
of the Lyapunov exponent, we choose the negative value 
for biological plausibility. Hence to lowest order we find 
the optimal phase resetting curve is Type II: 



A o (0) = - 



V^sin^Trf 



Va + Abn 2 + 16c7r 4 



(19) 



The next order correction does not appreciably change 
this result. To obtain the a 2 term, we must solve (fT5|) 
subject to ([FT)) . By the Fredholm Alternative, a solu- 
tion to the inhomogeneous problem exists if and only if 
the right-hand side of ([T5[) . call it r(6), is orthogonal to 
the nullspace of J* . However, since J is self-adjoint we 
simply solve for the value of v\ \ such that 



sm{2ir9)r{9)d0 = 0, 



namely, v± i = 0. 
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FIG. 1: In the case where the second derivative is left unconstrained, the optimal PRC deviates from a pure cosine function as 
the noise amplitude a increases. Parameters are a=l, b=l, c=0. 



Imposing periodic boundary conditions on the result- 
ing equation yields the solution 



Ai(0) = Cisin(27r6 



V2tt sin(27T0) sin(47r6 



(a - lUcir 4 )Va + Abir 2 + 16ctt 4 



As before, we use the constraint (fTTj) to obtain C\ = 0. 
Hence to order a 2 the optimal phase resetting curve is 
given by 



A(9) = 



+ 



V / 2sin(27r6') 
Va + 46tt 2 + 16c7T 4 
a 2 y/2n sin(27r6>) sin(47r6 



2 (a - 144c7r 4 )v'a + 46tt 2 + 16ctt 4 



(20) 



VI. CONSTRAINT PARAMETERS 

Let us next explore the influence of the constraint pa- 
rameters a, b and c, which we will allow to take on the val- 
ues of or 1. Of the seven nontrivial combinations, one 
has no periodic solution at all and is thus inadmissible. 
Four parameter choices give rise to the same optimum 
already found in Eq.(|20|). and two parameter combina- 
tions do not produce a unique solution but instead yield 
a family of solutions ranging smoothly from Type I to 
Type II. In this case, we explicitly find the minimizer of 
A among the family of solutions. 

All of the cases can be analyzed by examining Eq. (|T5f , 
the characteristic equation of C(A) = 0. For example, 



the case a = c = and b — 1 can have no periodic 
solution, since the polynomial (1 — v\fl)y 2 = has no 
nontrivial roots. 

The four parameter combinations that lead to Eq. ((20|) 
are those in which a = 1. In these cases we have 



^i,o + (1 



bfi,o)y 2 + cvi^y 4 



0. 



If c 7^ 0, the polynomial is 4th degree having four dis- 
tinct roots; if c = the polynomial is quadratic with two 
distinct roots. In each case we can set y = 2iri and solve 
uniquely for as discussed above. 

The case c = (while a = 1) deserves further attention 
for another reason. In this regime, the optimal PRC 
becomes sensitive to the noise amplitude a as illustrated 
in Fig. |T]). To understand why the curve deforms, let us 
focus on the extrema of Eq. ([20|) . which are given by the 
zeros of the derivative: 



A'(0) 



2y/2n 



cos(2tt6>) 



Va + 4&7T 2 + 16C7T 4 
n- 2 7T / 

cos(47T#) sin(27r6 



1 



144C7T 4 



cos(2ir9) sin(47r6 



In this form we clearly see that the unperturbed extrema 
(when cr = 0) occur at 6 = 1/4 and 3/4, while deforma- 
tion due to noise is on the order of er 2 7r/(a — 144c7r 4 ). 
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FIG. 2: When the first derivative is unconstrained while the second derivative is constrained, Euler-Lagrange optimization 
produces a family of candidates for the minimizer of the Lyapunov exponent ranging smoothly from Type II to Type I as the 
parameter K ranges from to 1. For negative K (dashed), the curves do not represent biologically plausible PRCs. Parameters 
are a = 0, b = 1, c = 1. 



More specifically, when c ^ this quantity is O(er 2 10~ 4 ) 
so that the weak noise in our model (a <C 1) has negligi- 
ble effect. However when c = 0, this quantity is 0(a 2 ), 
so that even relatively small magnitude noise can have a 
noticible impact on the shape of the optimal PRC. 

Another interesting situation arises in the two cases 
where a = 0, c = 1 and b is arbitrary. Here the charac- 
teristic equation has a double root at y = 0: 



(1 - bv lfi )y 2 + v lfl y 4 



0. 



After accounting for the boundary conditions, we have a 
superposition of two independent solutions 

A o (0) = C 3 (l - cos(2tt0)) + C 4 sin(27r6»). 



The constraint ([16]) eliminates only one degree of free- 
dom, leaving a family of solutions as candidates for the 
optimum: 



A o (0) = K 



1 — cos(27r6 



-yJ\-K 2 



-Air 2 ) 

sin(27r0) 
^27r 2 (6 + 47r 2 7' 



(21) 



where the remaining degree of freedom K has been nor- 
malized to range between —1 and 1. See Fig. @. 

Combining Eq.JU) for the Lyapunov exponent with 
Eq.© for the steady-state phase distribution, we insert 



Eq. ({2"Tj) to obtain the following expression: 

1 cr 4 (AK 4 + 10K 2 + 1) 



5 + 47T 2 4 47r 2 (6 + 4vr 2 ) 3 ' 

where we have set a — 0, c = 1. Note that we needed to 
carry out the expansion of A to cr 4 in order to discover 
the dependence on K. 

Since the derivative of A with respect to K has only one 
real root at K = 0, where a minimum occurs, the Type 
II curve remains the optimal PRC even in this case. 



VII. 



NUMERICAL VERIFICATION 



We would like to independently verify the accuracy 
of the optimal PRC ([2U|) derived via perturbation expan- 
sion by numerically solving the Euler-Lagrange equations 
(fTTjl and (|12p with periodic boundary conditions. Unfor- 
tunately, the resulting system is singular and therefore 
very difficult to solve numerically. Instead we substi- 



tute the approximation P(9) 
Euler-Lagrange functional ([S] 



= l + %-A(0)A'(0) into the 
to obtain a new functional 



£ A"A(1 + ^A(0)A'(0) 



+v x [aA 2 + b{A') 2 + c(A") 2 - l] d6 = 0, 
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FIG. 3: The magnitude of the optimal PRC depends on the whether or not the second derivative is constrained. The numerical 
solution (open circles) and the analytic result (solid lines) coincide. Parameters are a = 1,6 = 1 and a — 0.05. 



which gives rise via Eq.© to the 4th order boundary 
value problem 



A< 4 > = 



-2A" - 2aAv 1 + 2bN'v x 



V 2 - 3AA'A"o- 2 



When c = 0, we similarly obtain a 2nd order boundary 
value problem. 

Using the numerical integration package XPPAUT, we 
are able to achieve excellent agreement with our analyti- 
cal approximation. Fig. ([3]) illustrates numerical and an- 
alytic solutions in the case where c = 1 and where c = 0. 
Note that imposing a constraint on the second derivative 
of A results in an optimal PRC of much smaller magni- 
tude. 

In Fig.ffJI we find good agreement between the analytic 
and numerical results even for the regime in which a = 1 , 
c = and PRC shape is sensitive to noise amplitude. 
The numerical simulation deforms with increasing a just 
as the analytic approximation does. 



Discussion 

In this paper we have used perturbation theory and 
the calculus of variations to analyze the rate at which 
neurons can synchronize when subjected to common in- 
puts. We treat the inputs as "noise," that is, as if they 
are delta-correlated with no structure. Real neuronal in- 
puts do have correlational structure, however, so that the 



expression for the rate of synchronization (the Lyapunov 
exponent) is more complex. Indeed, in previous work 
[15| we have shown that the temporal characteristics of 
the noise can also have an effect on how rapidly neurons 
synchronize. In that work, we asked the reverse ques- 
tion: given a particular PRC, what correlation time for 
the noise minimizes the Lyapunov exponent? 

Suppose that we use some signal that is not white noise 
but still has zero mean and is stationary. Then the phase 
satisfies 

| = 1 + A( W ) 
where £(t) is the input. The Lyapunov exponent is 



A := lim — 

T^oo T 



dt. 



By using an approximation of 0(t) as in |20 | we may 
be able to obtain a functional for A depending on 
and A, and from this apply similar methods to estimate 
the optimal shape of the PRC given the statistics of the 
inputs. 

Optimization has been applied to other aspects of neu- 
ral oscillators. Moehlis, et al. [2l[ asked the following 
question. Consider the scalar oscillator model: 



dO 



f(9)+A(e)I(t). 



(Note that if f(9) — 1, we have Eq.©, the case consid- 
ered in this paper.) Suppose the neuron fired at t = and 
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FIG. 4: When the second derivative is unconstrained, the optimal PRC shape deforms with increasing noise. The numerical 
solution (open circles) and the analytic result (solid lines) are in good agreement. Parameters are a = 1, 6 = 1, c = 0. 



we desire it to fire again at time T > 0. What is the min- 
imum stimulus, I(t) (which, say, minimizes I(t) 2 dt) 
to do this? Moehlis, et al. [21| write the Euler-Lagrange 
equations for this optimization problem and then assume 
that I(t) is small in order to use perturbation methods. 
A related issue is the "optimal stimulus" [13] for produc- 
ing a spike in a neuron, and for neural oscillators this has 
been answered in 12311. 



APPENDIX: AN EXISTENCE PROOF 

On the interval [0,1], the phase resetting curve A is 
necessarily at the endpoints and possibly at interior 
points as well. As a result, we have a singular differential 
equation for the steady state distribution of phases P, 
derived earlier as Eq.© and repeated here: 



J= -P+ — A(AP)'. 



(A.l) 



However we will now see that Eq. (|A.l[) does indeed have 
a solution despite the singularities. 

Suppose A(9) 7^ in the open interval (a,b) C [0, 1], 
while A(a) = A(6) = 0. In this way, we will be able to 
apply our proof to the entire domain [0, 1] in a piecewise 
fashion; for example, if A(x) = sin(27ra;), then a = and 
b = 1/2, or a — 1/2 and 6=1. In the following we 
will assume, without loss of generality, that A(8) > in 
(a, 6). 



Let us begin by rewriting the differential equation as 
an integral equation. Define Q{x) := A(x)P(x). Then 
Eq. (|A.l[) becomes 



Q' 



2Q 



-2J 



a 2 A 2 a 2 A' 
We now introduce an integrating factor; let 



(A.2) 



z(x) :-- 



ds 



A 2 ( S )' 



where c € (o, b) is fixed. Observe that, as x approaches 
a from above we eventually have x < c, and hence z{x) 
approaches +oo. Likewise, as x approaches b from below, 
z(x) approaches — oo. 
Eq. (IA.2j) now becomes 

(e z ^Q)' = 
Integrating both sides gives 



2J 



,z(x) 



-z{x) 



K 



A(t) 



dt 



(A.3) 



where if is a constant of integration that will be deter- 
mined below. 

We see from Eq.(jAT]) that P(a) = P{b) = J. 
Therefore a solution exists iff lim x ^ a + Q(x) / A(x) = 
lim x ^ b - Q(x)/ A(x) = J. Let us first consider the right 
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endpoint and assume for now that the limit 

e*(t) 



lim 

x — >6~ 

exists. Let us compute 

,. Q{x) 2J 
i?-A(x) = 

2.7 



A(t) 



dt = L 



(A.4) 



lim 



AW 



K 



e z(t) 



dt 



lim 



AW 



A^e 2 ^) 



and note that when we set K ~ L, both numerator and 
denominator tend to as x — > b~. Thus we can use 
L'Hopital's rule to obtain 



lim 

x — >b~ 



Q(x) 



2. J 



lim 



/A(x) 



A(x) a 2 x^b- A(x)z'(x)e z(x '> + A'(x)e z ( x ) 



2J 



lim 



-1/A(x) 



^ Hm 
2 J / - 2 

7^" 



-^s^A(a;)+A'(a;) 
-1 



+ A(x)A'(x) 



(A.5) 



Now let us return to the assumption we made and ob- 
serve that the integral in Eq. (|A.4p is not improper af- 
ter all. Rewriting the integrand of (|A.4|) such that both 
numerator and denominator go to infinity, we can use 
L'Hopital's rule again to see that the integrand goes to 
zero: 



The last equality follows since A' is bounded and 
0. Hence our assumption was justified. 



lim. 



*b- e 



Now let us rewrite Eq. (|A.3[) , incorporating our knowl- 
edge from Eq. (|A.4p . namely that K = L: 



Q(x) = 




x e z (*) 



c A(i) 



dt 



It remains to show that lim x _ >a + Q(x)/ A(x) = J. We 
will prepare to use L'Hopital's rule once again by writing 



Q(x) 

Aa+ A(x) 



lim 



2J ^ e~*W 
a 2 x™+ A(x) J x A(t) 



dt 



2J .. I 
lim 



x A(t) 



dt 



, x (A.6) 

a 2 x^a+ A(x)e z W 

Since e z ^ tends to infinity as x approaches a from above, 
by L'Hopital's rule the denominator of (|A.6j) also tends 
to infinity: 



,z(x) 



x-,a+ 1/A(x) 



2 , e z W/A(x) 2 
-— hm 



a 2 x ^a+ A'(x)/A(x) 2 

2 e z ^1 
— hm — — — — 

a 2 x^a+ A'(x) 



The numerator of Eq. (|A.6[) tends to infinity as well since 



e z(t) 

lim -ttt 

t-6- A(t) 



1/A(t) 



lim 

t^b- e" z (*) 



-awaw! 

t-i- e- 2 (*)/A(t) 2 



lim -A'(t)e 
0. 



b e z ^ 
x A(t) 



A > 



71/ 



when M = max{A(x) : x G [0, 1]}, and the latter integral 
is clearly unbounded as x approaches a. Therefore we can 
apply to (|A.6[) a similar calculation to that in (|A.5P and 
conclude that \im x ^ a + Q(x)/A(x) = J as desired. 
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